/* contrib/earthdistance/earthdistance.c */

#include "postgres.h"

#include <math.h>

#include "utils/geo_decls.h"	/* 对于点 */

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

PG_MODULE_MAGIC;

/* 地球的半径以法定英里为单位。 */
static const double EARTH_RADIUS = 3958.747716;
static const double TWO_PI = 2.0 * M_PI;


/******************************************************
 *
 * degtorad - 将度转换为弧度
 *
 * arg: double, 以度为单位的角度
 *
 * returns: double, 以弧度表示的相同角度
 ******************************************************/

static double fc_degtorad(double fc_degrees)
{
	return (fc_degrees / 360.0) * TWO_PI;
}

/******************************************************
 *
 * geo_distance_internal - 点之间的距离
 *
 * args:
 *	 一对点 - 对于每个点，
 *	   x坐标是格林威治以西的经度（单位：度）
 *	   y坐标是赤道以北的纬度（单位：度）
 *
 * returns: double
 *	 两点之间在地球表面上的距离（单位：英里）
 ******************************************************/

static double fc_geo_distance_internal(Point *fc_pt1, Point *fc_pt2)
{
	double		fc_long1,
				fc_lat1,
				fc_long2,
				fc_lat2;
	double		fc_longdiff;
	double		fc_sino;

	/* 将度转换为弧度 */

	fc_long1 = fc_degtorad(fc_pt1->x);
	fc_lat1 = fc_degtorad(fc_pt1->y);

	fc_long2 = fc_degtorad(fc_pt2->x);
	fc_lat2 = fc_degtorad(fc_pt2->y);

	/* 计算经度差 - 想要 < 180 度 */
	fc_longdiff = fabs(fc_long1 - fc_long2);
	if (fc_longdiff > M_PI)
		fc_longdiff = TWO_PI - fc_longdiff;

	fc_sino = sqrt(sin(fabs(fc_lat1 - fc_lat2) / 2.) * sin(fabs(fc_lat1 - fc_lat2) / 2.) +
				cos(fc_lat1) * cos(fc_lat2) * sin(fc_longdiff / 2.) * sin(fc_longdiff / 2.));
	if (fc_sino > 1.)
		fc_sino = 1.;

	return 2. * EARTH_RADIUS * asin(fc_sino);
}


/******************************************************
 *
 * geo_distance - 点之间的距离
 *
 * args:
 *	 一对点 - 对于每个点，
 *	   x坐标是格林威治以西的经度（单位：度）
 *	   y坐标是赤道以北的纬度（单位：度）
 *
 * returns: float8
 *	 两点之间在地球表面上的距离（单位：英里）
 ******************************************************/

PG_FUNCTION_INFO_V1(geo_distance);

Datum geo_distance(PG_FUNCTION_ARGS)
{
	Point	   *fc_pt1 = PG_GETARG_POINT_P(0);
	Point	   *fc_pt2 = PG_GETARG_POINT_P(1);
	float8		fc_result;

	fc_result = fc_geo_distance_internal(fc_pt1, fc_pt2);
	PG_RETURN_FLOAT8(fc_result);
}
